library(EBImage)
library(ggplot2)
library(cowplot)
library(magrittr)


setwd("~/Michael_Ratz/cell_segmentation/")

Detection of EGFP expressing cells

This is probably not the most accurate segmentation algorithm, but hopefully it does the job for this particular application. I’ll try to go through every step of the process and maybe we can modify it later if something looks fishy.

read image data

First, let’s just read the olig2 tif and crop out a smaller region that we can use for segmentation. We also apply a normalization step to adjust the intensity values so that the highest values are close to 1 (see ?normalize).

olig2 <- "data/olig2.tif"

cells <- readImage(olig2)
cells <- normalize(cells)
cells_crop = cells[2000:3000, 2000:3000]
display(cells_crop, method = "raster")

Apply 2D convolution filter

Obviously we can see the cells already, but we have quite a lot of background that could interfere with the segmentation. We can apply a 2D convolution filter (see ?filter2) that can be used as a low-pass filter to remove noise from our image. First, we need to create a filter that will be used to compute the fast 2D FFT convolution product.

f = makeBrush(9, shape = 'disc', step = FALSE)

# Divide by sum
f = f/sum(f)

cells_crop_filtered <- filter2(cells_crop, filter = f)
display(cells_crop_filtered, method = "raster")

# Define filter function for later
filter_cells <- function(im, brush.size = 9, verbose = FALSE) {
  if (verbose) cat("Applying 2D convolution filter to image ... \n")
  f = makeBrush(brush.size, shape = 'disc', step = FALSE)
  f = f/sum(f)
  imfiltered <- filter2(im, filter = f)
  return(imfiltered)
}

intensity correction + segmentation

The changes are very subtle at this resolution, but you can see that the image has been slightly blurred.

The next step is to use the filtered version for thresholding. I came up with my own strategy here and I don’t know how applicable this method will be to different stainings. Anyways, below is a short explanation of the algorithm:

filter: run a 2D convolution filter on original image
correction: subract filtered image from original image (this should hopefully highlight the fluorescence signal) and normalize intensity values
threshold: calculate a threshold that is 2 standard deviations above the mean instensity and then segment out pixels above this threshold
cleaning: remove small object based on shape area
watershed: separate objects that are merged

What you typically would do is to segment based on intensities that are higher in the original image than in the filtered image + some threshold. This method usually works well, but I’ve noticed that you need to define a custom threshold for each data type. Instead, this method seems to be easier to adjust automatically to fit the data at hand.

We have already computed step 1 so the next code block will show how to run step 2 and 3 and later we’ll go though step 4 and 5. The histogram shows the distribution of instensity values for the corrected image and the dotted line shows the threshold used for segmentation.

# Correct image
cells_crop_corrected <- (cells_crop - cells_crop_filtered) %>% normalize()

# Define correction function for later
correct_cells <- function(im, imfiltered, verbose = FALSE) {
  if (class(im) != "Image") stop(paste0("Invalid input format of im: ", class(im)))
  if (class(imfiltered) != "Image") stop(paste0("Invalid input format of imfiltered: ", class(imfiltered)))
  if (verbose) cat("Correcting image ... \n")
  imcorrected <- (im - imfiltered) %>% normalize()
  return(imcorrected)
}

# Threshold image
thr <- mean(cells_crop_corrected) + 2*sd(cells_crop_corrected)
hist(cells_crop_corrected)
abline(v = thr, lty = "longdash")

cells_crop_th <- cells_crop_corrected > thr

# Define threshold function for later
threshold_cells <- function(imcorrected, nsd = 2, verbose = FALSE) {
  if (class(imcorrected) != "Image") stop(paste0("Invalid input format: ", class(imcorrected)))
  if (verbose) cat("Thresholding cells ... \n")
  thr <- mean(imcorrected) + nsd*sd(imcorrected)
  imthreshold <- imcorrected > thr
  return(imthreshold)
}

Below is a comparison between the original image, the corrected image and the segmented image. Note that the autofluorescent background is removed in the corrected image, e.g. the array spots.

par(mfrow = c(3, 1), mar = c(0, 0, 0, 0))

# Plot original image
display(cells_crop, method = "raster")
text(x = 10, y = 20, label = "original", adj = c(0,1), col = "orange", cex = 1.5)

# Plot corrected image
display(cells_crop_corrected, method = "raster")
text(x = 10, y = 20, label = "corrected", adj = c(0,1), col = "orange", cex = 1.5)

# Plot thresholded image
display(cells_crop_th, method = "raster")
text(x = 10, y = 20, label = "segmented", adj = c(0,1), col = "orange", cex = 1.5)

This method seems to work quite well at least for the this cropped window. We can also overlay the result to outline the cells in the original image.

cells_crop_colored <- rgbImage(red = cells_crop, green = cells_crop, blue = cells_crop)
cells_crop_th <- bwlabel(cells_crop_th)
cells_crop_colored <- paintObjects(cells_crop_th, cells_crop_colored, col = "#FFA500")
display(cells_crop_colored, method = "raster")
text(x = 10, y = 20, label = "outlined cells [Olig2]", adj = c(0,1), col = "orange", cex = 1.5)

Neun data

This method is not 100% perfect but hopefully it’ll be decent enough for it’s intended purpose. I would suspect that it fails in areas where the cells overlap reach other more. Lets’ see what it looks like if we use the Neun image instead.

Neun <- "data/neun.tif"

cells <- readImage(Neun)
cells <- normalize(cells)
cells_crop = cells[2000:3000, 2000:3000]
cells_crop_filtered <- filter2(cells_crop, filter = f)
cells_crop_corrected <- (cells_crop - cells_crop_filtered) %>% normalize()
thr <- mean(cells_crop_corrected) + 2*sd(cells_crop_corrected)
cells_crop_th <- cells_crop_corrected > thr
cells_crop_colored <- rgbImage(red = cells_crop, green = cells_crop, blue = cells_crop)
cells_crop_th <- bwlabel(cells_crop_th)
cells_crop_colored <- paintObjects(cells_crop_th, cells_crop_colored, col = "#FFA500")
display(cells_crop_colored, method = "raster")
text(x = 10, y = 20, label = "outlined cells [Neun]", adj = c(0,1), col = "orange", cex = 1.5)

Compute features

Now that we have defined a segmented image we can extract features from it. The function computeFeatures allows us to extract features such as total area (in pixels) and x, y coordinates among other things. The shape computation takes a lot of time to run so we might want to split the image into smaller batches to take down the computational time.

fts.shape <- computeFeatures.shape(x = cells_crop_th)
fts.moment <- computeFeatures.moment(x = cells_crop_th)

Here’s an explanation of the shape features: - s.area: area size (in pixels) - s.perimeter: perimeter (in pixels) - s.radius.mean: mean radius (in pixels) - s.radius.sd: standard deviation of the mean radius (in pixels) - s.radius.max: max radius (in pixels) - s.radius.min: min radius (in pixels)

And the momnet features: - m.cx: center of mass x (in pixels) - m.cy: center of mass y (in pixels) - m.majoraxis: elliptical fit major axis (in pixels) - m.eccentricity: elliptical eccentricity defined by sqrt(1-minoraxis^2/majoraxis2). Circle eccentricity is 0 and straight line eccentricity is 1. - m.theta: object angle (in radians)

We can explore some of these featuresif we want to, but what would be nice is if we can set a threshold to remove cells that are too small. I’m not sure what would be a good way of doing this but let’s explore the options.

gg <- as.data.frame(fts.shape)
ggm <- reshape2::melt(gg)

## No id variables; using all as measure variables

ggplot(data = ggm, aes(variable, value, fill = variable)) +
  geom_violin() +
  geom_jitter(size = 0.3, alpha = 0.5) +
  facet_wrap(~variable, scales = "free", ncol = 6) +
  ggtitle("cell shape features [Neun]")

Now we can of course also use the coordinates to plot the cell positions using other plot function s.a. ggplot.

gg <- as.data.frame(cbind(fts.shape, fts.moment))

ggplot(data = gg, aes(m.cx, 1000 - m.cy)) + 
  geom_point(size = 0.5) +
  ggtitle("cell distribution")

Let’s see what happens if we use the parameter features to filter out small cells. Now if you look closely you will see that the very smallest cells have dissapeared (we can also use the cell area).

par(mfrow=c(3, 1))
inds <- which(fts.shape[, "s.perimeter"] <= 5)
cells_crop_th_clean <- rmObjects(x = cells_crop_th, index = inds)
cells_crop_clean_colored <- rgbImage(red = cells_crop, green = cells_crop, blue = cells_crop)
cells_crop_th <- bwlabel(cells_crop_th)
cells_crop_clean_colored <- paintObjects(cells_crop_th_clean, cells_crop_clean_colored, col = "#FFA500")
display(cells_crop_colored, method = "raster")
text(x = 10, y = 20, label = "outlined cells [Neun]", adj = c(0,1), col = "orange", cex = 1.5)
display(cells_crop_clean_colored, method = "raster")
text(x = 10, y = 20, label = "outlined cells cleaned [Neun]", adj = c(0,1), col = "orange", cex = 1.5)
display(cells_crop, method = "raster")
text(x = 10, y = 20, label = "original image", adj = c(0,1), col = "orange", cex = 1.5)

# Define cleaning function for later
clean_cells <- function(imthreshold, ftr = "s.area", thr = 5, do.fast = FALSE, verbose = FALSE) {
  if (class(imthreshold) != "Image") stop(paste0("Invalid input format: ", class(imthreshold)))
  if (verbose) cat("Cleaning up unwanted speckles ... \n")
  imthreshold <- bwlabel(imthreshold)
  if (do.fast & ftr == "s.area") {
    inds <- which(table(imthreshold)[-1] < thr)
  } else if (do.fast & ftr != "s.area") (
    stop(paste0("Invalid option ", ftr, " for ftr. Only ftr = 's.area' can be used for fast shape extraction."))
  ) else {
    fts.shape <- shape_extractor(x = imthreshold)
    stopifnot(ftr %in% colnames(fts.shape))
    inds <- which(fts.shape[, ftr] <= thr)
  }
  imclean <- rmObjects(x = imthreshold, index = inds)
  return(imclean)
}

Split merged objects

The last step (5) will be to split objects that are merged. For this purpose we can use the watershed algorithm. Here’s the explanation of watershed from the EBImage package:

The algorithm identifies and separates objects that stand out of the background (zero). It inverts the image and uses water to fill the resulting valleys (pixels with high intensity in the source image) until another object or background is met. The deepest valleys become indexed first, starting from 1.

The plot below shows all the segmented cells and the color indicate uniqueness. If you look closely you will see that some longer shapes have now been split into multiple shapes by the watershed algorithm.

cells_crop_th_watershed <- watershed(x = EBImage::distmap(cells_crop_th), tolerance = 0.1)
display(colorLabels(cells_crop_th_watershed), method = "raster")
text(x = 10, y = 20, label = "labelled cells [Neun]", adj = c(0,1), col = "orange", cex = 1.5)

# Define watershed function for later
watershed_cells <- function(imclean, tol = 0.1, verbose = FALSE) {
  if (class(imclean) != "Image") stop(paste0("Invalid input format: ", class(imclean)))
  if (verbose) cat("Applying watershed ... \n")
  imwatershed <- watershed(x = EBImage::distmap(imclean), tolerance = tol)
  return(imwatershed)
}

Now we have a pretty decent workflow for segmenting cells so let’s combine the 5 steps into one function.

SegmentCells <- function(
  impath,
  crop.window = NULL,
  brush.size = 9,
  nsd = 3,
  feature = "s.area",
  feature.threshold = 5,
  do.fast = FALSE,
  tolerance = 0.1,
  return.all = FALSE,
  verbose = FALSE
) {
  if (!file.exists(impath)) stop(paste0("File ", impath, " does not exist \n"))
  cells <- readImage(impath)
  cells <- EBImage::normalize(cells)
  
  if (!is.null(crop.window)) {
    if (!length(crop.window) == 4 & class(crop.window) %in% c("numeric", "integer")) stop("Invalid crop window \n")
    cells <- cells[crop.window[1]:crop.window[2], crop.window[3]:crop.window[4]]
  }
  
  cells_filtered <- filter_cells(cells, brush.size, verbose) # 1. filter 
  cells_corrected <- correct_cells(cells, cells_filtered, verbose) # 2. correct 
  cells_th <- threshold_cells(cells_corrected, nsd, verbose) # 3. threshold
  cells_clean <- clean_cells(cells_th, feature, feature.threshold, do.fast, verbose) # 4. clean
  cells_split <- watershed_cells(cells_clean, tolerance, verbose) # 5. watershed
  
  if (return.all) {
    return(list(cells, cells_filtered, cells_corrected, cells_clean, cells_split))
  } else {
    return(cells_split)
  }
}

Apply segmentation

Now that we have a complete segmentation algoprithm we can apply it to the three different scans.

img.files <- c("data/olig2.tif", "data/egfp.tif", "data/neun.tif")

segmented <- list()
for (i in 1:3) {
  print(img.files[i])
  segmented[[i]] <- SegmentCells(
    impath = img.files[i], 
    crop.window = c(500, 4300, 230, 3900), 
    nsd = 2, 
    feature.threshold = 4, 
    do.fast = TRUE, 
    verbose = TRUE)
}

Find overlapping signals

The next step will be to find out what wignals are overlapping, but first we can just visualize the results from the three segmented images. Since we have three targets, we can encode them directly in RGB channels. You can already see some overlapping cells. We can ignore the cells along the edges for now, these will not be present when we run the whole image at once.

cat("Total number of oligodendrocytes: ", length(table(segmented[[1]])), "\n")

## Total number of oligodendrocytes:  13958

cat("Total number of EGFP: ", length(table(segmented[[2]])), "\n")

## Total number of EGFP:  3473

cat("Total number of neurons: ", length(table(segmented[[3]])), "\n")

## Total number of neurons:  37190

par(mfrow = c(3, 1))
display(1 - segmented[[1]], method = "raster")
text(x = 10, y = 20, label = "Olig2", adj = c(0, 1), col = "orange", cex = 1.5)

display(1 - segmented[[2]], method = "raster")
text(x = 10, y = 20, label = "EGFP", adj = c(0, 1), col = "orange", cex = 1.5)

display(1 - segmented[[3]], method = "raster")
text(x = 10, y = 20, label = "Neun", adj = c(0, 1), col = "orange", cex = 1.5)

im <- rgbImage(red = segmented[[1]], green = segmented[[2]], blue = segmented[[3]])
display(im, method = "raster")
text(x = 10, y = 20, label = "combined stainings", adj = c(0, 1), col = "orange", cex = 1.5)

With these segmented images it should be pretty straightforward to extract cells with overlapping signals.

neuron + EGFP = pixels with positive signal for both Neun and EGFP
oligodendrocyte + EGFP = pixels with positive signal for both Olig2 and EGFP
unknown + EGFP = pixels with positive signal for EGFP but negative signal for Olig2 and Neun

oligodendrocyte <- EBImage::channel(im, "red") > 0 & EBImage::channel(im, "green") > 0
neuron <- EBImage::channel(im, "blue") > 0 & EBImage::channel(im, "green") > 0
unknown <- EBImage::channel(im, "green") > 0 & EBImage::channel(im, "blue") == 0 & EBImage::channel(im, "red") == 0
olig2_neun <- EBImage::channel(im, "red") > 0 & EBImage::channel(im, "blue")

par(mfrow = c(5, 1))
display(dilate(oligodendrocyte), method = "raster")
text(x = 10, y = 20, label = "oligodendrocyte + EGFP", adj = c(0,1), col = "orange", cex = 1.5)
display(dilate(neuron), method = "raster")
text(x = 10, y = 20, label = "neurons + EGFP", adj = c(0,1), col = "orange", cex = 1.5)
display(dilate(unknown), method = "raster")
text(x = 10, y = 20, label = "unknown + EGFP", adj = c(0,1), col = "orange", cex = 1.5)
oligneu <- rgbImage(red = dilate(oligodendrocyte), green = dilate(neuron), blue = dilate(unknown))
display(oligneu, method = "raster")
text(x = 10, y = 20, label = "(unknown [blue],  olgiodendrocytes[red] \n and neurons [green]) + EGFP", adj = c(0,1), col = "orange", cex = 1.5)
display(dilate(olig2_neun), method = "raster")
text(x = 10, y = 20, label = "overlapping Neun and Olig2", adj = c(0,1), col = "orange", cex = 1.5)

We will have to relabel the new shapes to get an estimate on the number cells of the three categories.

par(mfrow = c(1, 3))
fts.shape.neurons <- computeFeatures.shape(bwlabel(neuron))
cat("Total number estimated of Egfp+ neurons: ", nrow(fts.shape.neurons), "\n")

## Total number estimated of Egfp+ neurons:  1635

hist(fts.shape.neurons[, 1], breaks = 20)

fts.shape.oligodenrocytes <- computeFeatures.shape(bwlabel(oligodendrocyte))
cat("Total number estimated of Egfp+ oligodendrocytes: ", nrow(fts.shape.oligodenrocytes), "\n")

## Total number estimated of Egfp+ oligodendrocytes:  505

hist(fts.shape.oligodenrocytes[, 1], breaks = 20)

fts.shape.unknown <- computeFeatures.shape(bwlabel(unknown))
cat("Total number estimated of Egfp+ unknown: ", nrow(fts.shape.unknown), "\n")

## Total number estimated of Egfp+ unknown:  5843

hist(fts.shape.unknown[, 1], breaks = 20)

If we sum up these values, the total number of overlapping shapes (cells) is larger than the total number of detected nuclei. I noticed that when taking the overlapping set of pixels, you will often have partial overlaps. 1 Egfp cell can for example overlap with two different Olig2 cells which will generate two overlapping shapes while in reality you have only 1 or 0 overlapping signals. For this reason I think we have to be quite conservative and come up with some way to determine what is an overlapping signal and what is not.

In the image below I have plotted the Olig2 and Egfp signals and even though most of the time you see a quite decent overlap, there are cases where overlap looks pretty bad (yellow color = intersecting cells).

par(mfrow = c(1, 3))
display(im[400:700, 400:700, 1:2], method = "raster")
text(x = 10, y = 20, label = "Olig2 [red] Egfp [green]", adj = c(0,1), col = "orange", cex = 2)
display(EBImage::channel(im, "green")[400:700, 400:700], method = "raster")
text(x = 10, y = 20, label = "Egfp", adj = c(0,1), col = "orange", cex = 2)
display(EBImage::channel(im, "red")[400:700, 400:700], method = "raster")
text(x = 10, y = 20, label = "Olig2", adj = c(0,1), col = "orange", cex = 2)

imt <- EBImage::channel(im, "red")[400:700, 400:700]

I have written a small function to deal with this problem and select cells based on an overlap criteria instead. This is to make sure that shapes that are split into multiple new shapes are not all included in the output. Only the shape with the highest overlap is kept and the overlap also have to be above a specified threshold (e.g. 50%). The overlap between two shapes is defined as the area of the intersect divided by the area of the smallest shape.

estimate intersect of overlapping cells
for cells with multiple overlaps, keep only the top hit
for each pair of overlapping shapes A and B, estimate the overlap as intersect(A, B)/min(A, B)
return cells with an overlap of at least 50%

ima <- EBImage::channel(im, "red")[400:700, 400:700]
imb <- EBImage::channel(im, "green")[400:700, 400:700]

# Define overlap function
overlap_fkn <- function(x, y, i) {
  return(i/pmin(x, y))
}

OverlapImages <- function(
  ima, 
  imb, 
  overlap.min = 0.5, 
  return.merged = FALSE, 
  return.indices = FALSE, 
  verbose = FALSE
) {
  
  if (length(dim(ima)) == 3 | length(dim(imb)) == 3) {
    stop("Invalid dims for ima or imb")
  }
  ima <- bwlabel(ima); imb <- bwlabel(imb)
  # Summarize overlap across images
  if (verbose) cat(paste0("Calculating intersect between images \n"))
  d <- data.frame(a = as.numeric(ima), b = as.numeric(imb))
  d <- table(d) %>% as.matrix()
  d <- d[-1, -1]
  
  # Extract area for image a and image b
  area_a <- table(ima)[-1]
  aa <- setNames(as.numeric(area_a), names(area_a))
  area_b <- table(imb)[-1]
  ab <- setNames(as.numeric(area_b), names(area_b))
  if (verbose) cat(paste0("Finished calculating intersect between ", length(area_a), " shapes in image a and ", length(area_b), " shapes in image b \n"))
  
  # Collect indices for overlaping shapes
  intersect_ab <- which(d > 0, arr.ind = T)
  indices <- setNames(data.frame(apply(do.call(rbind, (lapply(1:nrow(intersect_ab), function(i) {
    inds <- intersect_ab[i, ]
    c(rownames(d)[inds[1]], colnames(d)[inds[2]], d[inds[1], inds[2]])
  }))), 2, function(x) as.character(x)), stringsAsFactors = F), nm = c("ind.a", "ind.b", "intersect"))
  
  # Collect shape areas, indices and intersect
  df <- data.frame(area.a = aa[indices$ind.a], area.b = ab[indices$ind.b], 
                   inda = as.integer(indices$ind.a), indb = as.integer(indices$ind.b), 
                   intersect = as.numeric(indices$intersect), stringsAsFactors = F)
  
  # Calculate overlap
  if (verbose) cat(paste0("Calculating overlap between images using shape intersect \n"))
  df$overlap <- overlap_fkn(x = df$area.a, y = df$area.b, i = df$intersect)
  df$keep <- df$overlap > overlap.min
  df <- subset(df, keep)
  
  # Clean up image a nd image b
  ima_clean <- rmObjects(ima, index = setdiff(as.integer(names(table(ima)[-1])), df$inda))
  imb_clean <- rmObjects(imb, index = setdiff(as.integer(names(table(imb)[-1])), df$indb))
  
  # Should the clean images be merged?
  if (return.merged) {
    imres <- ima_clean | imb_clean
  } else {
    imres <- ima_clean & imb_clean
  }
  
  # return extra data
  if (return.indices) {
    return(list(bwlabel(imres), df))
  } else {
    return(bwlabel(imres))
  }
}

Now we can apply this overlap detection tool to pairs of segmented images. Remember that the segmented images are stored as different color channels in the object “im”, where red = Olig2, green = Egfp and blue = Neun.

The unknown Egfp cells are here defined as the Egfp cells with neither Olig2 or Neun overlap. This still means that there could be some overlap, just smaller than 50%.

library(zeallot)

c(oligodendrocyte, Olig2_inds) %<-% OverlapImages(EBImage::channel(im, "red"), EBImage::channel(im, "green"), return.indices = T)
c(neuron, Neun_inds) %<-% OverlapImages(EBImage::channel(im, "blue"), EBImage::channel(im, "green"), return.indices = T)
c(olig2_neun, Olig2_Neun_inds) %<-% OverlapImages(EBImage::channel(im, "red"), EBImage::channel(im, "blue"), return.indices = T)

unknown <- rmObjects(EBImage::channel(im, "green"), index = as.numeric(c(Olig2_inds$indb, Neun_inds$indb)))

par(mfrow = c(4, 1))
display(dilate(oligodendrocyte), method = "raster")
text(x = 10, y = 20, label = "oligodendrocyte + EGFP", adj = c(0,1), col = "orange", cex = 1.5)
display(dilate(neuron), method = "raster")
text(x = 10, y = 20, label = "neurons + EGFP", adj = c(0,1), col = "orange", cex = 1.5)
display(dilate(unknown), method = "raster")
text(x = 10, y = 20, label = "unknown + EGFP", adj = c(0,1), col = "orange", cex = 1.5)
display(dilate(olig2_neun), method = "raster")
text(x = 10, y = 20, label = "overlapping Neun and Olig2", adj = c(0,1), col = "orange", cex = 1.5)

fts.shape.neurons <- computeFeatures.shape(neuron)
fts.moment.neurons <- computeFeatures.moment(neuron)
cat("Total number of estimated Egfp+ neurons: ", nrow(fts.shape.neurons), "\n")

## Total number of estimated Egfp+ neurons:  1328

fts.shape.oligodenrocytes <- computeFeatures.shape(oligodendrocyte)
fts.moment.oligodenrocytes <- computeFeatures.moment(oligodendrocyte)
cat("Total number of estimated Egfp+ oligodendrocytes: ", nrow(fts.shape.oligodenrocytes), "\n")

## Total number of estimated Egfp+ oligodendrocytes:  362

fts.shape.olig2_neun <- computeFeatures.shape(olig2_neun)
fts.moment.olig2_neun <- computeFeatures.moment(olig2_neun)
cat("Total number of overlapping Neun-Olig2: ", nrow(fts.shape.olig2_neun), "\n")

## Total number of overlapping Neun-Olig2:  1040

fts.shape.unknown <- computeFeatures.shape(unknown)
fts.moment.unknown <- computeFeatures.moment(unknown)
cat("Total number of unknown: ", nrow(fts.shape.unknown), "\n")

## Total number of unknown:  1853

cat("Total number of Egfp+ cells after cleaning: ", sum(nrow(fts.shape.neurons) + nrow(fts.shape.oligodenrocytes) + nrow(fts.shape.unknown)))

## Total number of Egfp+ cells after cleaning:  3543

There are still a couple of cells too many (3543 - 3473 = 70). These are most likely overlapping Olig2, Neun and Egfp which have been counted twice. By removing the overlapping Olig2 a nd Neun we can probably get rid of these, but these are quite few so maybe it’s good enough?

We can now extract the feature coordinates and labels to have a data.frame format that is a bit easier to work with

df <- rbind(setNames(cbind(data.frame(fts.moment.neurons[, 1:2]), "neuron"), nm = c("x", "y", "celltype")),
            setNames(cbind(data.frame(fts.moment.oligodenrocytes[, 1:2]), "oligodendrocyte"), nm = c("x", "y", "celltype")),
            setNames(cbind(data.frame(fts.moment.unknown[, 1:2]), "unknown"), nm = c("x", "y", "celltype")))

ggplot(df, aes(x, 3000 - y, color = celltype)) + 
  geom_point(size = 0.5) +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = "black"), 
        legend.text = element_text(colour = "white")) +
  scale_color_manual(values = c("neuron" = "green", "oligodendrocyte" = "red", "unknown" = "blue")) +
  facet_wrap(~celltype, ncol = 1)

pre-segmented data

And here are the pre-segmented images. First, let’s load them an invert the intensity so that cells have a value of 1. One thing I noticed with this segmentation is that there are artificial “large cells” in regions of high cell density. But maybe that doesn’t make such a huge difference.

I also noticed that the images are not binarized, so I’ll have to threshold them. I set the intensity threshold to 0.5.

seg.files <- c("data/binarized_aligned/olig2-06.png", 
               "data/binarized_aligned/neun-05.png",
               "data/binarized_aligned/egfp-04.png")

Olig2 <- 1 - readImage(seg.files[1])[, , 1]
colorMode(Olig2) <- "Grayscale"
Olig2 <- Olig2 > 0.5
Olig2 <- bwlabel(Olig2)

Egfp <- 1 - readImage(seg.files[3])[, , 1]
colorMode(Egfp) <- "Grayscale"
Egfp <- Egfp > 0.5
Egfp <- bwlabel(Egfp)

Neun <- 1 - readImage(seg.files[2])[, , 1]
colorMode(Neun) <- "Grayscale"
Neun <- Neun > 0.5
Neun <- bwlabel(Neun)

im <- rgbImage(red = Olig2, green = Egfp, blue = Neun)
display(im, method = "raster")
text(x = 10, y = 20, label = "combined stainings", adj = c(0, 1), col = "orange", cex = 1.5)

Now we can run the overlap function to estimate how many cells are labelled as Neuron or Oligodndrocyte.

c(oligodendrocyte, Olig2_inds) %<-% OverlapImages(Olig2, Egfp, return.indices = T)
c(neuron, Neun_inds) %<-% OverlapImages(Neun, Egfp, return.indices = T)
c(olig2_neun, Olig2_Neun_inds) %<-% OverlapImages(Neun, Olig2, return.indices = T)

unknown <- rmObjects(Egfp, index = c(Olig2_inds$indb, Neun_inds$indb))

par(mfrow = c(4, 1))
display(dilate(oligodendrocyte), method = "raster")
text(x = 10, y = 20, label = "oligodendrocyte + EGFP", adj = c(0, 1), col = "orange", cex = 1.5)
display(dilate(neuron), method = "raster")
text(x = 10, y = 20, label = "neurons + EGFP", adj = c(0,1), col = "orange", cex = 1.5)
display(dilate(unknown), method = "raster")
text(x = 10, y = 20, label = "unknown + EGFP", adj = c(0,1), col = "orange", cex = 1.5)
display(dilate(olig2_neun), method = "raster")
text(x = 10, y = 20, label = "overlapping Neun and Olig2", adj = c(0, 1), col = "orange", cex = 1.5)

fts.shape.Egfp <- computeFeatures.shape(Egfp)
fts.moment.Egfp <- computeFeatures.moment(Egfp)
cat("Total number of Egfp+ cells: ", nrow(fts.shape.Egfp), "\n")

## Total number of Egfp+ cells:  2912

fts.shape.neurons <- computeFeatures.shape(neuron)
fts.moment.neurons <- computeFeatures.moment(neuron)
cat("Total number of estimated Egfp+ neurons: ", nrow(fts.shape.neurons), "\n")

## Total number of estimated Egfp+ neurons:  879

fts.shape.oligodenrocytes <- computeFeatures.shape(oligodendrocyte)
fts.moment.oligodenrocytes <- computeFeatures.moment(oligodendrocyte)
cat("Total number of estimated Egfp+ oligodendrocytes: ", nrow(fts.shape.oligodenrocytes), "\n")

## Total number of estimated Egfp+ oligodendrocytes:  146

fts.shape.olig2_neun <- computeFeatures.shape(olig2_neun)
fts.moment.olig2_neun <- computeFeatures.moment(olig2_neun)
cat("Total number of overlapping Neun-Olig2: ", nrow(fts.shape.olig2_neun), "\n")

## Total number of overlapping Neun-Olig2:  197

fts.shape.unknown <- computeFeatures.shape(unknown)
fts.moment.unknown <- computeFeatures.moment(unknown)
cat("Total number of unknown: ", nrow(fts.shape.unknown), "\n")

## Total number of unknown:  1937

cat("Total number of Egfp+ cells after cleaning: ", sum(nrow(fts.shape.neurons) + nrow(fts.shape.oligodenrocytes) + nrow(fts.shape.unknown)))

## Total number of Egfp+ cells after cleaning:  2962

Date

date()

## [1] "Wed Mar 25 13:47:08 2020"

Session Info

sessionInfo()

## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] zeallot_0.1.0  magrittr_1.5   cowplot_1.0.0  ggplot2_3.2.1 
## [5] EBImage_4.26.0
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.2          plyr_1.8.4          compiler_3.6.1     
##  [4] pillar_1.4.2        bitops_1.0-6        tools_3.6.1        
##  [7] digest_0.6.22       evaluate_0.14       tibble_2.1.3       
## [10] gtable_0.3.0        lattice_0.20-38     pkgconfig_2.0.3    
## [13] png_0.1-7           rlang_0.4.1         yaml_2.2.0         
## [16] parallel_3.6.1      xfun_0.10           withr_2.1.2        
## [19] stringr_1.4.0       dplyr_0.8.3         knitr_1.25         
## [22] fftwtools_0.9-8     htmlwidgets_1.5.1   locfit_1.5-9.1     
## [25] grid_3.6.1          tidyselect_0.2.5    glue_1.3.1         
## [28] R6_2.4.0            jpeg_0.1-8          rmarkdown_1.16     
## [31] reshape2_1.4.3      purrr_0.3.2         scales_1.0.0       
## [34] htmltools_0.4.0     BiocGenerics_0.30.0 abind_1.4-5        
## [37] assertthat_0.2.1    colorspace_1.4-1    tiff_0.1-5         
## [40] labeling_0.3        stringi_1.4.3       RCurl_1.95-4.12    
## [43] lazyeval_0.2.2      munsell_0.5.0       crayon_1.3.4

Cell segmentation

Ludvig Larsson

3/22/2020